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Abstract 

A  steady  two-dimensional  computational  model  for  a  proton  exchange  membrane  (PEM)  fuel  cell  is  presented.  The  model  accounts  for 
species  transport,  electrochemical  kinetics,  energy  transport,  current  distribution,  and  water  uptake  and  release  in  the  catalyst  layer.  The 
governing  differential  equations  are  solved  over  a  single  computational  domain,  which  consists  of  a  gas  channel,  gas  diffusion  layer,  and 
catalyst  layer  for  both  the  anode  and  cathode  sides  of  the  cell  as  well  as  the  solid  polymer  membrane.  The  model  for  the  catalyst  regions  is 
based  on  an  agglomerate  geometry,  which  requires  water  species  to  exist  in  both  dissolved  and  gaseous  forms  simultaneously.  Data  related  to 
catalyst  morphology,  which  was  required  by  the  model,  was  obtained  via  a  microscopic  analysis  of  a  commercially  available  membrane 
electrode  assembly  (MEA).  The  coupled  set  of  differential  equations  is  solved  with  the  commercial  computational  fluid  dynamics  (CFD) 
solver,  CFDesign™,  and  is  readily  adaptable  with  respect  to  geometry  and  material  property  definitions.  The  results  show  that  fuel  cell 
performance  is  highly  dependent  on  catalyst  structure,  specifically  the  relative  volume  fractions  of  gas  pores  and  polymer  membrane 
contained  within  the  active  region  as  well  as  the  geometry  of  the  individual  agglomerates. 
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1.  Introduction 

In  recent  years,  the  proton  exchange  membrane  (PEM) 
fuel  cell  has  received  a  great  deal  of  attention,  from  the 
automotive  industry  in  particular,  as  a  candidate  for  near¬ 
future  power  generation  applications.  PEM  fuel  cells  are 
particularly  suited  to  automotive  applications  primarily 
because  of  their  relatively  low  operating  temperature,  effi¬ 
ciency,  and  high  power  density  that  is  comparable  to  existing 
internal  combustion  technology.  However,  like  many  emer¬ 
ging  technologies,  PEM  fuel  cells  must  overcome  certain 
engineering  and  economic  obstacles  if  they  are  ever  to 
become  commercially  viable.  In  short,  PEM  fuel  cells  must 
become  more  efficient  and  lower  in  cost.  Improvements  in 
cell  design  and  materials  development  can  help  to  achieve 
these  goals.  There  are  two  approaches  to  improving  cell 
design  and  materials  development.  The  first  is  to  design  and 
build  test  cells  and  evaluate  their  performance.  This  method 
can  yield  useful  information,  but  it  is  also  costly  and  time 
consuming.  In  addition,  although  it  is  relatively  easy  to 
evaluate  a  cell’s  total  current  and  voltage  and  from  these 
generate  a  polarization  curve,  it  is  much  more  difficult  to 
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evaluate  operating  parameters  in  situ,  which  is  vital  to 
understanding  how  a  design  performs  over  a  range  of 
operating  conditions.  In  order  to  obtain  data  such  as  reactant 
concentration  profiles,  membrane  hydration,  and  tempera¬ 
ture  distributions  inside  the  cell,  it  is  more  convenient  to 
simulate  the  fuel  cell  with  a  mathematical  model. 

In  the  past,  most  mathematical  models  focused  on  the 
cathode  side  of  the  fuel  cell  only;  the  reason  being  that  the 
cathode  activation  overpotential  is  the  single  largest  source  of 
inefficiency  in  the  fuel  cell.  These  models  also  typically 
include  the  membrane,  as  it  contributes  to  ohmic  overpoten¬ 
tial,  and  the  gas  diffusers.  The  catalyst  layer  itself  has 
generally  been  modeled  as  an  interface  and  denoted  the  point 
at  which  source  terms  for  species  consumption  or  production 
were  applied.  The  work  by  Springer  et  al.  [1]  laid  the 
foundation  for  many  future  numerical  models.  In  their  paper, 
they  present  a  one-dimensional  model,  the  principal  focus  of 
which  is  water  transport  through  the  membrane.  Their  model 
shows  that  the  net  water  transport  per  proton  is  much  less  than 
the  measured  electro-osmotic  drag  coefficient  for  a  fully 
hydrated  membrane,  which  indicates  the  presence  of  other 
important  water  transport  mechanisms.  Bernardi  and  Ver- 
brugge  [2]  present  a  model  with  the  focus  on  species  trans¬ 
port,  electrochemistry,  and  catalyst  utilization  in  which  they 
conclude,  among  other  things,  that  the  transport  of  gases 
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dissolved  in  the  polymer  phase  of  the  electrodes  affects 
catalyst  utilization  and  limits  cell  performance.  This  model 
is  based  on  a  pseudo-homogeneous  catalyst  layer  structure  in 
which  the  layer  has  no  pores  through  which  gas  can  be 
transported  to  the  reaction  sites.  Rather,  the  reactants  must 
dissolve  into  the  polymer  phase  and  diffuse  through  the  non- 
porous  catalyst  layer  to  reach  the  reaction  sites.  This  approach 
can  be  computationally  advantageous  as  the  equations 
describing  transport  in  the  catalyst  layer  are  relatively  simple. 
The  disadvantage  is  that  the  pseudo-homogenous  catalyst 
layer  structure  does  not  allow  for  transport  of  reactants  in  the 
gas  phase  within  this  layer.  Research  has  shown  that  the 
catalyst  layer  is  porous  and  that  reactants  can  be  transported 
throughout  in  the  gas  phase  [3,4].  The  model  of  Broka  and 
Ekdunge  [4]  represents  one  of  the  first  applications  of  an 
agglomerate  catalyst  layer  structure  to  a  PEM  model.  In  this 
model,  Broka  and  Ekdunge  show,  through  microscopic  ana¬ 
lysis,  that  the  catalyst  layer  is  made  up  of  clumps  of  carbon- 
supported  Pt  catalyst  surrounded  by  a  thin  layer  of  Nation  and 
separated  by  pores.  These  clumps  are  referred  to  as  agglom¬ 
erates.  The  principal  difference  between  this  type  of  model 
and  the  pseudo-homogenous  model  is  that  in  the  agglomerate 
model,  reactants  can  move  in  the  gas  phase  rather  than  solely 
as  a  dissolved  species  through  the  catalyst  layer.  Broka  and 
Ekdunge  also  show  that  the  agglomerate  model  is  better 
suited  to  modeling  fuel  cell  behavior  at  high  current  densities, 
where  concentration  overpotential  becomes  dominant.  In 
their  model,  they  assume  that  the  gaseous  reactant  concen¬ 
tration  is  constant  across  the  catalyst  layer.  In  the  model 
presented  in  this  paper,  the  reactant  concentration  is  treated  as 
variable. 

In  recent  years,  modeling  efforts,  both  our  own  [5-7]  and 
that  of  others  (e.g.  [8-12]),  have  increased  in  complexity. 
Most  current  models  are  multi-dimensional,  include  mass 
transport  of  multiple  species,  and  include  the  entire  fuel  cell 
geometry  from  one  gas  channel  to  the  other.  Some  models 
also  account  for  two-phase  flow,  which  is  necessary  when 
modeling  the  effects  of  liquid  water  production  and  flooding. 
Zhou  and  Liu  [8]  developed  a  three-dimensional  fuel  cell 
model  based  on  the  earlier  work  of  Gurau  et  al.  [9].  As  part 
of  the  work  they  show  how  the  model  solution  is  affected 
when  transitioning  from  a  two-dimensional  to  a  three- 
dimensional  geometry.  Their  model  focuses  on  species 
transport  as  well  as  current  and  temperature  distribution. 


To  solve  the  model,  they  partitioned  the  solution  domain  into 
three  coupled  regions,  which  is  a  more  involved  process  for 
the  user  than  is  the  single  domain  approach  used  in  this 
work.  In  addition,  water  uptake  and  release  within  the 
polymer  portion  of  the  catalyst  layers  is  neglected.  Um 
et  al.  [10]  present  a  transient,  three-dimensional  model 
and  show  that  an  interdigitated  flow  field  can  help  to  reduce 
mass  transfer  limitations.  They  use  a  single  domain  solution 
approach  and  neglect  water  uptake  and  release  in  the  catalyst 
layer.  Shimpalee  et  al.  [11]  also  present  a  three-dimensional 
model  but  do  not  include  transport  through  the  catalyst  layer 
as  it  is  modeled  as  an  interface.  They  show  that  the  direction 
of  water  transport  through  the  membrane  can  affect  current 
density  distribution  patterns.  A  model  including  two-phase 
flow  is  presented  by  Natarajan  et  al.  [12].  They  show  that 
liquid  water  buildup  in  the  cathode  has  a  substantial  influ¬ 
ence  on  cell  performance.  Their  model  is  two-dimensional 
and  focuses  on  liquid  water  transport  in  the  porous  media  of 
the  gas  diffuser.  The  catalyst  layer  is  an  interface  and  only 
the  cathode  side  of  the  cell  is  included.  It  should  be  noted 
that  with  the  exception  of  Broka  and  Ekdunge,  all  of  the 
models  mentioned  above  either  treat  the  catalyst  layer  as  a 
pseudo-homogeneous  film  or  simply  as  an  interface. 

This  paper  describes  a  mathematical  model  that  simulates 
the  transport  of  gaseous  species,  energy,  protonic  current, 
and  water  dissolved  in  the  polymer  phase  of  the  catalyst 
layers  and  membrane.  The  catalyst  layer  is  modeled  as 
having  an  agglomerate  structure,  and  the  effect  of  catalyst 
layer  structure  on  cell  performance  is  also  examined.  The 
model  is  based  on  the  finite  element  method  and  is  for¬ 
mulated  as  a  single  solution  domain  problem.  Consequently, 
interface  conditions  between  the  individual  fuel  cell  ele¬ 
ments  need  not  be  specified,  and  the  only  boundary  condi¬ 
tions  required  are  at  the  outer  surfaces  of  the  model. 

2.  Model  development 

Fig.  1  shows  the  solution  domain  of  the  model.  In  the 
anode  and  cathode  gas  channels,  fuel  and  oxidant  flow  along 
the  surface  of  the  membrane  electrode  assembly.  In  these 
regions,  the  flow  is  considered  to  be  laminar.  Reactants  move 
from  the  gas  channels  into  the  gas  diffusion  layers  (GDL) 
which  consist  of  a  thin  sheet  of  carbon  paper,  the  purpose  of 
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Fig.  1.  Solution  domain. 
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which  is  to  evenly  distribute  the  reactants  across  the  catalyst 
layers  and  provide  an  electrical  connection  with  the  collector 
plate  (not  shown).  The  catalyst  layers  are  among  the  more 
important  parts  of  the  fuel  cell  as  it  is  within  these  regions 
where  all  of  the  electrochemical  reactions  take  place.  A  large 
portion  of  the  work  presented  in  this  model  involves  the 
characterization  and  modeling  of  the  catalyst  layers.  The 
polymer  membrane,  which  is  assumed  to  be  impermeable  to 
reactant  gases,  transports  only  protons  and  dissolved  water. 
Both  ionic  conductivity  and  ionomer  content  in  the  MEA 
have  a  significant  of  impact  on  cell  performance. 

2.1.  Catalyst  layer  structure 

To  develop  an  accurate  description  of  the  catalyst  layer,  a 
5  cm2  MEA  was  purchased  from  ElectroChem  Inc.  Cross 
sections  of  the  MEA  were  prepared  and  analyzed  with  a 
scanning  electron  microscope  (SEM)  and  a  transmission 
electron  microscope  (TEM).  The  cross  sections  for  SEM 
analysis  were  prepared  by  three  methods  in  order  to  evaluate 
which  worked  best.  The  first  method  involved  freezing  the 
sample  in  liquid  nitrogen  and  fracturing  it.  This  technique 
provided  a  clean  section  for  viewing  but  also  caused  the  gas 
diffusion  layer  to  de-laminate  from  the  assembly.  The 
second  method  was  to  simply  cut  the  MEA  at  room  tem¬ 
perature  with  a  razor.  The  blade  was  pushed  directly  down 
on  top  of  the  MEA  as  opposed  to  cutting  across  it.  This  also 
resulted  in  a  clean  section,  but  there  was  some  dragging  of 
the  polymer  through  the  sample  as  it  is  fairly  malleable  at 
room  temperature.  For  the  third  method,  the  sample  was 
frozen  in  liquid  nitrogen  and  then  cut  as  per  the  second 
method.  This  provided  the  best  results  of  all. 

Preparation  of  the  cross  sections  for  TEM  analysis 
involved  encasing  the  sample  in  epoxy  and  then  using  a 
diamond  knife,  or  microtome,  to  slice  off  very  thin  sections 
(^70  nm)  for  viewing.  The  samples  were  cut  at  room 
temperature. 

The  structure  of  the  fuel  cell  MEA  is  shown  in  Fig.  2.  The 
region  labeled  A  in  this  image  is  the  gas  diffusion  layer. 
Region  B  is  the  catalyst  layer.  In  this  region  there  is  no 
carbon  paper,  only  polymer  encapsulated  catalyst  sites 
called  agglomerates.  Region  C  is  the  polymer  membrane. 

Fig.  3  shows  a  higher  magnification  view  of  the  catalyst 
layer  labeled  B  in  Fig.  2.  Fig.  3  along  with  Fig.  2  can  be  used 
to  estimate  the  mean  agglomerate  size,  the  thickness  of  the 
active  layer,  and  the  void  fraction  of  the  active  layer.  The 
agglomerate  size  and  void  fraction  were  evaluated  through 
the  use  of  image  analysis  and  enhancement  software.  By 
enhancing  the  SEM  and  TEM  images,  it  is  possible  to 
highlight  certain  catalyst  characteristics  of  interest  and  per¬ 
form  quantitative  analyses.  Geometric  information  about  the 
catalyst  structure  can  be  used  in  conjunction  with  manu¬ 
facturer’s  data  for  catalyst  loading  to  estimate  the  specific 
catalyst  area.  Pores  can  be  seen  in  the  catalyst  layer  on  the 
order  of  1-10  pm  in  size.  These  are  referred  to  as  macro¬ 
pores. 


Fig.  2.  Freeze  cut  cross  section  of  ElectroChem  MEA.  Image 
magnification  is  200  x. 

Fig.  4  shows  a  schematic  of  the  agglomerate  catalyst  layer 
geometry  used  in  this  model.  Reactant  gas  flows  through  the 
catalyst  macro-pores  and  then  dissolves  into  and  diffuses 
through  the  agglomerate  to  the  carbon-supported  platinum 
reaction  sites. 


20Pm  EHT  =  5.00  kV  WD  =  6  mm 


| - 1  Photo  No.  =  3995  Mag=  500  X 

Fig.  3.  Room  temperature  cut  section  of  MEA  showing  the  catalyst  layer. 
Image  magnification  is  500  x. 


Fig.  4.  Agglomerate  catalyst  geometry. 
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Fig.  5.  TEM  image  of  an  agglomerate.  Magnification  is  18,400x. 


Figs.  5  and  6  show  TEM  images  of  an  agglomerate.  In 
these  images,  it  is  possible  to  see  the  individual  Pt  catalyst 
particles  as  well  as  the  carbon  support.  As  indicated  in  Fig.  5, 
the  light  gray  area  between  the  carbon  and  Pt  particles  is 
ionomer.  The  white  areas  are  holes,  not  to  be  confused  with 
pores.  These  holes  result  from  the  pullout  of  carbon-sup- 
ported  catalyst  during  sample  preparation.  Indeed,  nano¬ 
scale  pores  are  absent  from  the  TEM  images.  This  is  an 
important  observation  as  it  indicates  that  gas  transfer  in  the 
porous  catalyst  is  facilitated  primarily  by  macro-pores  on  the 
order  of  1-10  pm  in  size. 

Fig.  6  is  a  higher  magnification  image  in  which  the  Pt 
catalyst  particles  are  well  defined.  These  Pt  particles  are 
nominally  3  nm  in  diameter.  However,  as  a  consequence  of 
the  technique  used  to  deposit  the  catalyst,  the  individual 
particles  have  a  tendency  to  clump  together.  When  the 
catalyst  particles  form  clumps,  there  is  a  reduction  in  the 
amount  of  surface  area  available  for  reaction,  which  can  be 
viewed  as  a  decrease  in  catalyst  utilization. 


Fig.  6.  TEM  image  of  an  agglomerate  showing  Pt  catalyst  particles  and 
clumping.  Image  magnification  is  485,500x. 

The  physical  parameters  used  in  the  model  can  be  found  in 
Table  5.  These  parameters,  some  of  which  were  determined 
by  the  microscopic  analyses  discussed  earlier,  must  be  used 
in  conjunction  with  the  appropriate  mathematical  relations 
to  accurately  model  fuel  cell  performance. 

2.2.  Governing  equations  and  closure  relations 

Eq.  (1)  in  Table  1  reflects  conservation  of  mass  for  all  gas 
species.  A  source  term  reflects  changes  in  the  overall  gas- 
phase  mass  due  to  consumption  or  production  of  gas  species 
via  reaction  and  mass  transfer  between  the  water  in  the  gas 
phase  and  that  dissolved  in  the  polymer.  The  total  density  of 
the  gas  mixture  is  expressed  as  the  sum  of  the  individual 
species  concentrations  multiplied  by  their  respective  molar 
masses. 

The  momentum  equations,  Eqs.  (2a)  and  (2b),  are 
expressed  as  the  Navier-Stokes  equations  in  vector  form, 
modified  with  a  source  term  to  account  for  Darcy  flow  in  the 


Table  1 


Governing  equations 


Governing  equation 

Vector  form  of  equation 

Equation 

Global  continuity 

V  ■  ( p8U )  =  5h2  +  So2  +  Syap  +  Sdiss 

i 

Momentum 

U  ■  V[psux }  =  -dP/dx  +  /l V2ux  +  SDarx 

2a 

U  ■  =  -dP/dy  +  /(V2uv  +  Soary 

2b 

Hydrogen  transport 

V  ■  (Dh2 Vch2)  —  SVch2  +  Sh2  =  0 

3 

Water  vapor  transport 

V  ■  (-D»H2oVc,,h2o)  —  »Vc„h2o  +  Svap  +  Sdiss  =  0 

4 

Oxygen  transport 

V  ■  (Dq2  Vco2)  —  mVco2  +  >So2  =  0 

5 

Nitrogen  transport 

V  ■  (£>n2Vcn2)  -  SVcn2  =  0 

6 

Dissolved  water  transport 

V  •  (DmVcm)  +  Sdrag  —  Sdiss  —  0 

7 

Potential  membrane 

V  ■  (a mV(/)m)  =  Sem 

8 

Thermal  energy 

V  •  [keffVT]  —  pCeffUS/T  +  iSom  +  Srev  +  *^act  —  Spl  —  Spc  —  ‘S'evap  =  0 

9 
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Table  2 
Source  terms 


Source  term 

Equation 

Application 

Darcy  flow  in  ^-direction,  Sbaix  — 

C/iux 

10 

Anode  and  cathode:  GDL,  catalyst 

Darcy  flow  in  y-direction,  Soary  = 

C/uiv 

11 

Anode  and  cathode;  GDL,  catalyst 

Hydrogen  consumption,  5h2  —  — (1/2F)BV 

12 

Anode  catalyst 

Water  vapor  production,  5Vap  —  (1/2F)  BV 

13 

Cathode  catalyst 

Vapor/dissolved  water  mass  transfer,  Sdiss  - 

^mass(^m  GH20) 

14 

Cathode  catalyst 

Oxygen  consumption,  5o2  =  —  (1/4F)BV 

15 

Cathode  catalyst 

Electro-osmotic  drag,  Sdrag  —  —(2.5/22 F)  (82/ dx)o(dcp / dx) 

16 

Anode  and  cathode:  catalyst,  membrane 

Protonic  current,  Sem  —  —BY 

17 

Anode  and  cathode:  catalyst 

Ohmic  heating,  Som  =  (^m9<l>m/dx)2(l/am 

) 

18 

Anode  and  cathode:  catalyst 

Reversible  heat,  Siev  —  F(BV/F) 

J2s°/n 

19 

Anode  and  cathode:  catalyst 

Activation  loss,  Sact  —  (</>c  —  (j)m ) 

BV 

20 

Anode  and  cathode:  catalyst 

Water  vaporization,  Spc  —  (BV/2F)/ifg 

21 

Cathode  catalyst 

Vapor/dissolved  phase  change,  Sevap  —  ^fg^diss 

22 

Anode  and  cathode:  catalyst 

Collector  plate  heat  sink,  Spi  -  [(7  -  TpiJ/ffigdiAgdi)  +  (fcouAcoii))]  1/fgdi 

23 

Anode  and  cathode;  GDL 

porous  regions  of  the  model.  The  Darcy  source  is  active  in 
the  GDL  and  catalyst  layers  only;  the  inertial  and  viscous 
terms  are  neglected  in  these  regions. 

The  gas  species  equations  are  given  in  Eqs.  (3)-(6).  The 
diffusion  coefficients  are  based  on  a  simplification  of  the 
Stefan-Maxwell  equations  and  are  modified  by  a  porosity 
factor  [9].  Each  of  these  equations  has  an  advective  term 
equal  to  the  product  of  velocity  and  concentration  gradient. 
The  general  form  of  these  equations  would  also  include  the 
complement  to  this  term,  the  product  of  concentration  and 
velocity  gradient.  In  preliminary  simulations,  the  velocity 
gradient  was  found  to  be  negligible  for  most  of  the  solution 
domain,  with  the  exception  of  the  inlet  and  exit,  and  for  this 
reason  the  product  of  concentration  and  velocity  gradient  is 
omitted  from  the  gas  species  equations.  The  source  terms  for 
hydrogen  and  oxygen  species,  Eqs.  (12)  and  (15)  in  Table  2, 
account  for  consumption  via  reaction.  The  source  term  for 
water  vapor  accounts  for  production  of  water  at  the  cathode, 
Eq.  (13). 

Water  exists  in  dissolved  form  within  the  polymer  mem¬ 
brane  and  a  portion  of  the  polymer  phase  of  the  catalyst 
layer.  Dissolved  water  is  transported  through  the  polymer  by 
diffusion  and  electro-osmotic  drag  only.  A  convective  term 
would  have  to  be  added  if  the  gas  pressures  in  the  anode  and 
cathode  were  different;  for  this  work  the  anode  and  cathode 
pressures  are  within  3%  of  each  other  and  so  the  convective 
transport  of  water  through  the  MEA  is  neglected.  A  source 
term  is  needed  to  account  for  mass  transfer  between  the 
dissolved  and  vapor  phases  within  the  catalyst  layer.  This 
term,  Eq.  (14),  is  applied  to  both  the  dissolved  water  species 
equation  and  to  the  water  vapor  equation.  Eq.  (7)  is  the 
transport  equation  for  dissolved  water. 

The  transport  of  protons  in  the  polymer  portions  of  the 
fuel  cell  is  described  by  Eq.  (8).  The  source  term,  Eq.  (17), 
represents  the  production/consumption  of  protons  via  the 
electrochemical  reactions  in  the  catalyst  layers.  The  rate  of 
the  electrochemical  reaction  is  described  by  the  Butler- 


Volmer  relation,  Eq.  (24).  The  electrical  potential  is  assumed 
to  be  constant  over  each  electrode.  It  is  set  to  zero  on  the 
anode  side  and  to  the  difference  between  the  cell  voltage  and 
the  open  circuit  voltage  on  the  cathode  side. 

The  energy  equation  is  expressed  by  Eq.  (9)  and  contains 
sources  for  ohmic  heating  due  to  ionic  resistance  Eq.  (18), 
reversible  heat  Eq.  (19),  heat  produced  via  activation  losses 
Eq.  (20),  and  heat  exchange  involved  in  the  phase  change  of 
water.  There  are  two  terms  related  to  phase  change  energy 
transfer.  The  first  accounts  for  the  energy  needed  to  vaporize 
the  water  produced  via  reaction  Eq.  (21).  One  of  the  assump¬ 
tions  of  this  model  is  that  water  is  produced  in  the  liquid  phase 
and  instantly  vaporized,  provided  that  the  cathode  stream  is 
not  saturated.  The  energy  required  for  the  vaporization  must 
then  be  included.  The  other  term  accounts  for  water  moving 
between  the  dissolved  and  vapor  phases  Eq.  (22).  The 
enthalpy  of  the  water  dissolved  in  the  polymer  is  assumed 
to  be  the  same  as  the  enthalpy  of  liquid  water. 

In  this  model,  the  collector  plates  are  not  included  in  the 
solution  domain.  In  an  actual  fuel  cell,  the  shoulder  of  the 
plates  would  be  in  contact  with  the  MEA  and  provide  a  low- 
resistance  pathway  for  heat.  The  source  term  Spi  Eq.  (23) 
approximates  the  amount  of  heat  that  would  move  through  the 
collector  plates  had  they  been  part  of  the  solution  domain. 

Table  3  contains  the  closure  relations  needed  to  complete 
the  mathematical  model.  Eq.  (24)  is  the  Butler- Volmer 
equation  as  expressed  for  an  agglomerate  catalyst  geometry 
[4],  The  first  term  in  the  denominator  effectively  sets  the 
maximum  flux  of  reactant  that  can  pass  through  the  polymer 
layer  surrounding  an  agglomerate.  This  term  is  based  on  a 
dissolved  gas  concentration  at  the  gas/polymer  interface 
given  by  a  Henry’s  law  relation,  Eq.  (26).  The  reactant  flux 
limit  establishes  the  limiting  current  of  the  cell.  The  second 
term  in  the  denominator  includes  the  reaction  rate,  Eq.  (25), 
and  is  similar  to  that  used  by  Broka  and  Ekdunge.  with  the 
exception  that  in  this  model  the  concentration  of  reactant 
(either  oxygen  or  hydrogen)  is  assumed  to  vary  across  the 
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Table  3 

Closure  relations 


Butler-Volmer,  BV  (A/nm3) 

BV  =  nF  /  ((<5/ a)/ C*Dmll)  +  1/kE 

24 

Reaction  rate,  k  (mol/mm3  s) 

k  =  A,ii0,e/nF(c*/ckfe{y[e^~'l,^F/RT  -  e~(^~K)^F/RT^ 

25 

Dissolved  gas  concentration,  c*  (mol/mm3) 

c *  =  hifiTck 

26 

Effectiveness,  E 

E  =  tanh  (mL)/mL 

27 

Thiele’s  modulus,  mL 

mE  —  Ly^k/(^c*Dmg) 

28 

Gas  species  diffusion  coefficient,  Dk  (mm2/s) 

Dk  =  [((1  -  ckMWk)/p)/J2(cj/cmDjk)\  A5 

Dm  =  1.3E  —  4e2416(1/303'-(i/7’) 

29 

Diffusion  coefficient  for  dissolved  water,  Dm  (mm2/s) 

31 

Protonic  conductivity,  om  (Q  mm)-1 

am  =  (0.0051397  -  0.000326)  e(126S/(1/303)-(1/J'))(TgTm) 

32 

Membrane  water  content,  X  (mol  H20/mol  SO3-) 

7  =  cm/((MWm/pm)  -  0.0126 cm) 

34 

width  of  the  catalyst  layer.  Hence  the  term  c*,  which  is  the 
dissolved  reactant  concentration,  is  not  constant  The  effec¬ 
tiveness,  given  as  Eq.  (27),  is  a  measure  of  how  readily 
reactants  diffuse  through  the  agglomerates  [4],  An  effec¬ 
tiveness  of  1.0  would  indicate  that  reactants  diffusing 
through  the  agglomerates  encounter  no  resistance.  Eq.  (29) 
is  an  effective  diffusion  coefficient  which  allows  for  an 
approximate  solution  to  the  multi-component  species  equa¬ 
tions  without  having  to  solve  the  full  set  of  Stefan-Maxwell 
equations  [9].  Eqs.  (31)  and  (32)  are  the  diffusion  coeffi¬ 
cients  for  water  in  Nafion  and  the  ionic  conductivity  of 
Nafion,  respectively  [1].  Eq.  (34)  relates  membrane  water 
content  to  the  dissolved  water  concentration  while  account¬ 
ing  for  the  swelling  of  the  membrane  with  water  sorption. 

Table  4  contains  a  listing  of  the  boundary  conditions 
required  by  the  model.  Boundary  conditions  are  specified 
at  the  inlets  and  outlets  to  both  gas  channels  while  the 
remaining  external  surfaces  of  the  solution  domain  are 
assumed  to  have  zero  flux  conditions. 

Although  the  model  presented  does  not  include  transient 
effects,  initial  conditions  (ICs)  are  specified  for  the  species 


equations  as  well  as  the  dissolved  water  equation.  Specify¬ 
ing  ICs  helps  to  stabilize  the  numerical  solution  and  also 
decreases  the  amount  of  time  needed  to  reach  a  converged 
solution.  The  ICs  for  the  species  equation  are  the  same  as  the 
inlet  boundary  conditions.  The  dissolved  water  equation  IC 
is  set  in  a  manner  to  be  consistent  with  the  water  vapor 
equation  to  which  it  is  coupled  by  Eq.  (14). 

In  addition  to  the  governing  equations  and  closure  rela¬ 
tions,  certain  values  for  material  properties  and  other  phy¬ 
sical  parameters  are  needed  to  complete  the  model.  The 
numerical  values  for  the  physical  parameters  used  in  the 
model  are  given  in  Table  5. 

2.3.  Numerical  considerations 

The  outer  surfaces  of  the  gas  channels  shown  in  Fig.  1  are 
bounded  by  the  collector  plates  (not  shown)  and  are 
impermeable  to  gases.  As  part  of  the  single  domain  for¬ 
mulation,  each  governing  equation  is  solved  throughout  the 
entire  domain,  even  if  the  equation  is  not  physically  valid  in 
every  region.  This  is  accomplished  using  a  variety  of 


Table  4 

Boundary  conditions  for  base  case 


Governing  equation 

Anode  inlet 

Anode  exit 

Cathode  inlet 

Cathode  exit 

Comments 

Momentum 

216  mm/s  average 
velocity,  parabolic  profile 

Fixed  pressure 

422  mm/s  average 
velocity,  parabolic  profile 

Fixed  pressure 

- 

Hydrogen  transport 

Ch2  —  8.57E  —  8  mol/mm3 

V  ■  cH,  =  0 

o 

II 

s? 

o 

> 

V  ■  ch2  =  0 

Initial  concentration 

in  cathode  is 

set  to  zero 

Oxygen  transport 

<1 

f-5 

o 

II 

o 

<1 

r> 

O 

II 

O 

c'o2  —  1.84E  —  8  mol/mm3 

V  ■  c02  =  0 

Initial  concentration 

in  anode  is 

set  to  zero. 

Water  vapor  transport 

cvh2o  =  8.65E  —  9  mol/mm3 

V  ■  c„h2o  —  0 

V  ■  cvh2o  =  7.98E  —  9  mol/mm3 

V  ■  c„h2o  =  0 

- 

Nitrogen  transport 

o 

II 

z 

G 

> 

<1 

Z 

II 

o 

cn2  =  6.94E  —  8  mol/mm3 

V  ■  cN2  =  0 

Initial  concentration 

in  anode  is 

set  to  zero 

Dissolved  water  transport 

o 

II 

E 

O 

> 

V  ■  Cm  =  0 

V  ■  Cm  =  0 

V  ■  Cm  =  0 

Initial  concentration 
based  on  equilibrium 
with  water 

vapor 

Membrane  potential 

V  ■  <!)m  =  0 

v  -K  =  0 

v  •  * m  =  0 

V  ■  K  =  0 

- 

Energy 

348.0  K 

V-  T  =  0 

348.0  K 

V  T  =  0 

Constant  temperature 
BC’s  along  gas 
channels 
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Table  5 

Physical  parameters 


Gas  channel  length,  Lgc  (cm)  6.6 

Gas  channel  width,  W  (mm)  2.0 

Gas  channel  height,  H  (mm)  2.0 

Gas  diffusion  layer  thickness,  rgdl  (mm)  0.175 

Catalyst  layer  thickness,  tcat  (mm)  0.100 

Membrane  thickness,  tm  (mm)  0.127 

Collector  thickness,  rcol  (mm)  3.0 

Shoulder  width,  w  (mm)  1.0 

Reference  temperature,  Tref  (K)  298 

Reference  pressure,  Pref  (atm)  1.0 

Fixed  charge  concentration,  cH+  (mol/m* * 3)  1,200.0 

Gas  diffusion  layer  void  fraction,  ig  0.5 

Pt/carbon  volume  fraction  in  the  catalyst  layer,  iptc  0.19 

Catalyst  layer  void  fraction,  rcat  0.09 


Current  Density,  A/cmA2 

Fig.  7.  Numerical  and  experimental  performance  curves  for  test  cell. 


Table  6 

Physical  conditions  for  fuel  cell  test 


Ionomer  volume  fraction  in  the  catalyst  layer,  iion 

0.72 

Test  parameter 

Value 

Air  inlet  relative  humidity,  Rhc  (%) 

65 

Fuel  inlet  relative  humidity,  Rha  (%) 

60 

Anode  inlet  pressure  (psig) 

24.9 

Faraday’s  constant,  F  (C/mol) 

96,485 

Cathode  inlet  pressure  (psig) 

25.5 

Permeability  to  air  of  gas  diffusion 

5.68  x  104 

Cell  temperature  (K) 

348 

layer,  K  (mm-2) 

Anode  volume  flow  rate  (SLPM) 

2.5 

Cathode  viscosity,  /rair  (Pa  s) 

1.0  x  10~5 

Cathode  volume  flow  rate  (SLPM) 

50 

Anode  viscosity,  (Pa  s) 

2.0  x  10~5 

Anode  relative  humidity  (%) 

65 

Specific  porous  area  of  the  catalyst  layer. 

8,000 

Cathode  relative  humidity  (%) 

60 

Mean  agglomerate  diameter,  Dagg  (mm) 

Agglomerate  characteristic  length,  L  (mm) 
Agglomerate  diffusion  parameter,  S/a  (mm2) 
Diffusivity  of  oxygen  in  the  polymer,  Dm  o2  (mm/s) 
Diffusivity  of  hydrogen  in  the  polymer,  Dm  h2  (mm/s) 
Reference  anode  exchange  current  density,  i‘o,a 
(A/mm2) 

Reference  cathode  exchange  current  density,  z‘o,c 
(A/mm2) 

Anodic  transfer  coefficient,  ya 
Cathodic  transfer  coefficient,  ac 
Oxygen  reference  concentration,  Co2,ref  (mol/mm3) 
Hydrogen  reference  concentration,  fH2,ref  (mol/mm3) 
Convective  mass  transfer  coefficient,  ftraass  (s  ') 
Solubility  coefficient  for  the  cathode,  hAc  ( K  '  ) 


Solubility  coefficient  for  the  anode,  /td  a  (K 
Open  circuit  voltage,  Voc  (V) 


6.1E-3 
1.02  x  10~3 
2.42  x  10~3 
1.15  x  10~4 
2.29  x  10~4 
5.0  x  10~6 * * * 

6.0  x  10~10 


12 

3/2 

4.55  x  10~ 
2.17  x  10_ 
100 

4.09E-4 

1.82E-3 

1.0 


numerical  techniques  [9,13].  The  domain  is  divided  into 
35  x  147  elements.  Mapped  meshing  was  used  to  maintain  a 
sufficient  mesh  density  throughout  the  model  domain.  A 
sensitivity  analysis  was  conducted  by  doubling  the  number 

of  elements  in  the  mesh.  The  solution  changed  on  average  by 
less  than  0.33%  and  so  was  assumed  to  be  mesh  independent. 


3.  Results  and  discussion 

Fig.  7  shows  an  experimental  polarization  data  set  for  a 

50  cm2  ElectroChem  cell  operating  at  the  conditions  given 

in  Table  6.  Also  plotted  in  Fig.  7  is  a  curve  showing  the  cell 

polarization  predicted  using  the  model  operating  under  the 
same  conditions.  The  numerical  results  are  consistent  with 

the  data  set  over  the  entire  operating  range. 

In  this  model,  the  transport  of  reactants  through  the 

catalyst  layer  occurs  along  two  parallel  pathways.  The  first 


is  gas  transport  by  diffusion  and  advection  through  the 
macro-pores  within  the  catalyst  layer  (shown  in  Fig.  2). 
As  the  reactants  pass  through  the  pores,  they  dissolve  into 
the  catalyst  agglomerates  and  diffuse  inward  through  the 
surrounding  polymer  to  the  carbon-supported  platinum 
reaction  sites.  For  a  fixed  volume  of  carbon-supported 
catalyst  particles,  a  change  in  the  catalyst  layer  void  fraction 
requires  a  corresponding  change  in  the  fraction  of  ionomer 
in  the  catalyst  layer.  Fig.  8  shows  cell  performance  as  void 
fraction  is  changed  at  a  cell  voltage  of  0.5  V. 

The  optimal  catalyst  layer  void  fraction  for  the  test  MEA 
at  this  voltage  is  0.04  or  4%.  The  results  in  Fig.  8  show  that 
for  void  fractions  less  than  0.04  the  cell  performance  drops 
rapidly.  This  is  due  to  rapidly  increasing  concentration 
overpotential.  As  the  void  fraction  and  permeability  decrease, 
reactant  transport  by  diffusion  and  advection  within  the 
catalyst  layer  drops  sharply.  The  effect  is  exacerbated  by 
the  fact  that  the  catalyst  layers  on  the  test  MEA  are  relatively 


Cell  Voltage  =  0.5  V 


0  0.05  0.1  0.15  0.2  0.25  0.3  0.35 

Catalyst  Void  Fraction 

Fig.  8.  Cell  performance  variation  with  catalyst  layer  void  fraction. 
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Fig.  9.  (a)  Oxygen  concentration  change  across  the  cathode  catalyst  layer, 
(b)  Oxygen  concentration  change  down  the  cathode  flow  channel,  9%  void 
fraction.  Inlet  conditions  as  per  Table  4. 


thick,  100  pm.  Catalyst  layer  thickness  also  plays  a  role  in  the 
performance  drop  off  for  void  fractions  greater  than  0.04.  As 
the  void  fraction  increases,  the  ionomer  fraction  must 
decrease.  Thus,  while  an  increase  in  porosity  improves  mass 
flow  and  reduces  concentration  overpotential,  it  also  increases 
the  ohmic  overpotential.  Since  the  catalyst  layers  are  so  thick 
(relatively  long  ionic  conduction  pathway),  the  test  cell  is 
particularly  sensitive  to  changes  in  ionic  conductivity. 

Fig.  9a  shows  the  variation  of  oxygen  concentration 
through  the  fuel  cell  (perpendicular  to  the  flow  direction) 
as  the  void  fraction  is  varied.  Fig.  9b  shows  the  oxygen 
variation  down  the  channel  at  the  catalyst  layer/membrane 
interface. 

Fig.  9a  shows  that  the  overall  decrease  in  concentration  of 
oxygen  through  the  catalyst  is  not  a  linear  function  of  void 
fraction.  As  the  void  fraction  changes  from  0.04  to  0.02,  the 
concentration  drop  increases  by  roughly  15%.  As  the  void 
fraction  changes  from  0.02  to  0.01  the  concentration  drop 
increases  by  approximately  44%.  This  can  be  attributed  to 
the  effect  of  permeability  on  reactant  transport.  Fig.  9b 
shows  that  at  the  base  case  conditions,  the  oxygen  concen¬ 
tration  drops  by  5%  down  the  channel.  The  non-linear 
portion  of  the  curve  can  be  attributed  to  developing  flow 
effects  and  higher  local  reaction  rates. 


Fig.  10.  Cell  performance  variation  with  characteristic  length  of  the 
catalyst  agglomerates. 

Fig.  10  shows  the  effect  of  the  characteristic  length  of  the 
catalyst  agglomerates  on  fuel  cell  performance. 

The  results  in  this  figure  indicate  that  the  size  of  the 
catalyst  agglomerates  can  significantly  affect  cell  perfor¬ 
mance.  As  agglomerate  size  increases,  the  reactants  must 
diffuse  through  a  greater  distance  to  reach  the  catalyst  sites, 
thus  causing  an  increase  in  concentration  overpotential.  The 
smallest  possible  characteristic  length  is  about  0.01  pm, 
which  is  the  ratio  of  the  volume  to  surface  area  for  a  single 
catalyzed  carbon  support  particle  with  an  assumed  spherical 
geometry  (~50  nm  diameter).  One  way  to  reduce  the  char¬ 
acteristic  length  is  to  create  a  catalyst  layer  with  smaller 
agglomerates.  Another  way  is  to  control  the  shape  of  the 
agglomerates,  i.e.  use  a  structure  with  a  low  volume  to 
surface  area  ratio,  such  as  a  ribbon.  For  the  test  cell  used 
in  validating  the  model,  the  average  characteristic  agglom¬ 
erate  length  was  determined  from  the  SEM  analysis  to  be 
about  1  pm. 

4.  Conclusions 

A  steady  state,  two-dimensional  model  of  a  proton 
exchange  membrane  fuel  cell  has  been  developed  for  the 
purpose  of  predicting  fuel  cell  behavior  over  a  range  of 
operating  conditions  and  for  use  as  a  design  tool.  Work 
currently  underway  will  expand  this  model  to  three-dimen¬ 
sional  and  incorporate  liquid  water  transport. 

Results  from  our  two-dimensional  model  were  presented 
and  validated  with  experimental  data.  Electron  microscopy 
was  found  to  be  an  effective  method  for  determining  certain 
physical  parameters  used  in  the  model  such  as  catalyst  void 
fraction,  agglomerate  size,  and  catalyst  layer  thickness.  The 
void  fraction  of  the  catalyst  layer  was  found  to  significantly 
influence  cell  performance.  Model  results  show  that  a  void 
fraction  of  0.04  in  the  catalyst  layer  is  optimal  for  the 
ElectroChem  MEA  on  which  the  model  is  based.  In  addition 
increasing  the  agglomerate  characteristic  length  leads  to  a 
decrease  in  cell  performance  primarily  due  to  increased 
diffusive  resistance  to  reactant  flow.  Control  of  catalyst  layer 
structure  at  the  microscopic  level,  particularly  void  fraction 
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and  characteristic  agglomerate  length,  could  lead  to  better 
cell  performance  in  the  high  current  density  region  where 
concentration  overpotential  is  most  significant. 
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